---
title: "Replication Code"
format: html
---

```{r}
if (!require("pacman")) install.packages("pacman")
pacman::p_load(readr,dplyr,ggplot2,lme4,stringr,tidyr,forcats,sandwich,lmtest,flexmix,mlogit,broom,knitr,stringr,numDeriv,flextable,broom,stringr,car)
```

# Utilities

```{r}
floor_dec <- function(x,digit=1){
  floor(x*10^digit)/10^digit
}
ceiling_dec <- function(x,digit=1){
  ceiling(x*10^digit)/10^digit
}
monotonically_increasing <- function(v){
  mi = TRUE
  for(i in 1:length(v)){
    if(i<length(v)){
      if(!all(v[i]<=v[(i+1):length(v)])){
        mi = FALSE
        break
      }
    }
  }
  return(mi)
}

monotonically_decreasing <- function(v){
  mi = TRUE
  for(i in 1:length(v)){
    if(i<length(v)){
      if(!all(v[i]>=v[(i+1):length(v)])){
        mi = FALSE
        break
      }
    }
  }
  return(mi)
}

FMR_eva <- function(object,label){
  eva_matrix <- data.frame(k=object@k,logLik=NA,AIC=NA,BIC=NA)
  for(model in object@models){
    if(model@k0!=model@k){
      next
    }
    eva_matrix$logLik[eva_matrix$k==model@k] <- model@logLik
    eva_matrix$AIC[eva_matrix$k==model@k] <- AIC(model)
    eva_matrix$BIC[eva_matrix$k==model@k] <- BIC(model)
  }
  eva_matrix
}
```

# Simulations

The function of Equation (3):

```{r}
redist_per <- function(e_i,e_j,e_max,ml,mh,alpha,beta,theta_0,r_wage=c("h","l")){
  if(r_wage=="h"){
    m_j <- mh
  }else{
    m_j <- ml
  }
  theta <- theta_0+beta*(e_j/e_max)
  t<-(((mh*e_i)-theta^(1/(alpha-1))*m_j*e_j)/
  (1+(theta)^(1/(alpha-1))))/(mh*e_i)
  return(ifelse(t>=0,t,0))
} 
```

```{r}
e_max <- 10
ml <- 1
mh <- 2
alpha <- 0.8

```

```{r}
sim_data <- rbind(cbind(seq(0,10,0.01)*ml,
sapply(seq(0,10,0.01),redist_per,e_i=10,e_max=e_max,ml,mh,alpha,beta=0.9,theta_0=0.6,r_wage="l"),ml,"Highly responsive to effort"),
cbind(seq(0,10,0.01)*mh,
      sapply(seq(0,10,0.01),redist_per,e_i=10,e_max=e_max,ml,mh,alpha,beta=0.9,theta_0=0.6,r_wage="h"),mh,"Highly responsive to effort"),
cbind(seq(0,10,0.01)*ml,
      sapply(seq(0,10,0.01),redist_per,e_i=10,e_max=e_max,ml,mh,alpha,beta=0,theta_0=0.8,r_wage="l"),ml,"Less responsive to effort"),
cbind(seq(0,10,0.01)*mh,
      sapply(seq(0,10,0.01),redist_per,e_i=10,e_max=e_max,ml,mh,alpha,beta=0,theta_0=0.8,r_wage="h"),mh,"Less responsive to effort"),
cbind(seq(0,10,0.01)*ml,
      sapply(seq(0,10,0.01),redist_per,e_i=10,e_max=e_max,ml,mh,alpha,beta=0,theta_0=0.001,r_wage="l"),ml,"Self-interested"),
cbind(seq(0,10,0.01)*mh,
      sapply(seq(0,10,0.01),redist_per,e_i=10,e_max=e_max,ml,mh,alpha,beta=0,theta_0=0.001,r_wage="h"),mh,"Self-interested"))

sim_data<-data.frame(sim_data)
colnames(sim_data) <- c("income","redistribution","wage","cue")
sim_data$income<-as.numeric(sim_data$income)
sim_data$redistribution<-as.numeric(sim_data$redistribution)*100
sim_data$wage<-factor(sim_data$wage,labels=c("Low (1 unit)","High (2 units)"))
sim_data$cue<-factor(sim_data$cue,levels=c("Highly responsive to effort","Less responsive to effort","Self-interested"))
```

Figure 1:
```{r}
ggplot(data = sim_data, aes(x = income,y=redistribution))+ 
  geom_line(aes(colour=wage,linetype=wage),linewidth = 1)+
  scale_x_continuous(name="Recipient\'s income (tokens)",limits=c(0,20))+
  scale_y_continuous(name="% of income redistributed")+
  scale_colour_manual(name = "Recipient\'s wage",values=c("#F2C45F","#1A80BB"))+
  scale_linetype_manual(name = "Recipient\'s wage",values=c("longdash","dashed"))+
  facet_wrap(~cue)+
  theme_bw(base_size = 14)+
  theme(legend.position="bottom")
ggsave("./Figures/Figure 1.pdf",width=180,height=112,units = "mm")
```

# The full effort information condition

Import MTurk data:
```{r}
ob_A_mTurk <- read_csv("Observable_MTurk.csv")
```
Import student data:
```{r}
ob_A_lab <- read_csv("Observable_students.csv")
```

Transform to the long data format:
```{r}
ob_A_mTurk_long <- ob_A_mTurk %>%
  mutate(income_cue=if_else(yo_l_01>=yo_h_01 & yo_l_02>=yo_h_02 &yo_l_03>=yo_h_03 &yo_l_04>=yo_h_04 &yo_l_05>=yo_h_05 &yo_l_06>=yo_h_06 &yo_l_07>=yo_h_07 &yo_l_08>=yo_h_08 &yo_l_09>=yo_h_09 &yo_l_10>=yo_h_10,1,0),
         effort_cue=if_else(yo_l_02>=yo_h_01 & yo_l_04>=yo_h_02 & yo_l_06>=yo_h_03 & yo_l_08>=yo_h_04 & yo_l_10>=yo_h_05,1,0),
         increasing = if_else(monotonically_increasing(c(yo_l_00,yo_l_01,yo_l_02,yo_l_03,yo_l_04,yo_l_05,yo_l_06,yo_l_07,yo_l_08,yo_l_09,yo_l_10))&monotonically_increasing(c(yo_h_00,yo_h_01,yo_h_02,yo_h_03,yo_h_04,yo_h_05,yo_h_06,yo_h_07,yo_h_08,yo_h_09,yo_h_10)),1,0),
         decreasing = if_else(monotonically_decreasing(c(yo_l_00,yo_l_01,yo_l_02,yo_l_03,yo_l_04,yo_l_05,yo_l_06,yo_l_07,yo_l_08,yo_l_09,yo_l_10))&monotonically_decreasing(c(yo_h_00,yo_h_01,yo_h_02,yo_h_03,yo_h_04,yo_h_05,yo_h_06,yo_h_07,yo_h_08,yo_h_09,yo_h_10)),1,0),
         constant = if_else(increasing==decreasing & increasing!=0,1,0),
         increasing = increasing - constant,
         decreasing=decreasing-constant,
         all_zero = if_else(constant == 1 & yo_l_00==0,1,0),
         pattern = case_when(increasing==1~"increasing",decreasing==1~"decreasing",all_zero==1~"all_zero",TRUE ~ "other"))%>%
  pivot_longer(cols=yo_l_00:yo_h_10,names_to = "question",values_to ="transfer") %>%
  mutate(wage_j = if_else(str_sub(question,4,4)=="l",0.1,0.2),
         effort_j = as.numeric(str_sub(question,-2)),
         income_j = wage_j*effort_j,
         income_post_transfer = earning-transfer,
         transfer_per=(transfer/earning)*100,
                                belief = case_when(str_sub(question,-2)=="00" & str_sub(question,4,4)=="l" ~believe_income_00,
                                                   str_sub(question,-2)=="01" & str_sub(question,4,4)=="l" ~believe_income_01,
                                                   str_sub(question,-2)=="02" & str_sub(question,4,4)=="l" ~believe_income_02,
                                                   str_sub(question,-2)=="03" & str_sub(question,4,4)=="l" ~believe_income_03,
                                                   str_sub(question,-2)=="04" & str_sub(question,4,4)=="l" ~believe_income_04,
                                                   str_sub(question,-2)=="05" & str_sub(question,4,4)=="l" ~believe_income_05,
                                                   str_sub(question,-2)=="06" & str_sub(question,4,4)=="l" ~believe_income_06,
                                                   str_sub(question,-2)=="07" & str_sub(question,4,4)=="l" ~believe_income_07,
                                                   str_sub(question,-2)=="08" & str_sub(question,4,4)=="l" ~believe_income_08,
                                                   str_sub(question,-2)=="09" & str_sub(question,4,4)=="l" ~believe_income_09,
                                                   str_sub(question,-2)=="10" & str_sub(question,4,4)=="l" ~believe_income_10,
                                                   str_sub(question,-2)=="00" & str_sub(question,4,4)=="h" ~believe_income_00,
                                                   str_sub(question,-2)=="01" & str_sub(question,4,4)=="h" ~believe_income_02,
                                                   str_sub(question,-2)=="02" & str_sub(question,4,4)=="h" ~believe_income_04,
                                                   str_sub(question,-2)=="03" & str_sub(question,4,4)=="h" ~believe_income_06,
                                                   str_sub(question,-2)=="04" & str_sub(question,4,4)=="h" ~believe_income_08,
                                                   str_sub(question,-2)=="05" & str_sub(question,4,4)=="h" ~believe_income_10,
                                                   str_sub(question,-2)=="06" & str_sub(question,4,4)=="h" ~believe_income_12,
                                                   str_sub(question,-2)=="07" & str_sub(question,4,4)=="h" ~believe_income_14,
                                                   str_sub(question,-2)=="08" & str_sub(question,4,4)=="h" ~believe_income_16,
                                                   str_sub(question,-2)=="09" & str_sub(question,4,4)=="h" ~believe_income_18,
                                                   str_sub(question,-2)=="10" & str_sub(question,4,4)=="h" ~believe_income_20,),
         ep4ew=(earning+income_j)/(num_correct+effort_j)*num_correct,
         t_ep4ew=if_else(earning-ep4ew>0,earning-ep4ew,0),
         b_ep4ew=if_else(((floor_dec(t_ep4ew)>0|ceiling_dec(t_ep4ew)>0)&transfer>0)|((floor_dec(t_ep4ew)==0|ceiling_dec(t_ep4ew)==0)&transfer==0),TRUE,FALSE),
         equality=(earning+income_j)/2,
         t_equality=if_else(earning-equality>0,earning-equality,0),
         b_equality=if_else(((floor_dec(t_equality)>0|ceiling_dec(t_equality)>0)&transfer>0)|((floor_dec(t_equality)==0|ceiling_dec(t_equality)==0)&transfer==0),TRUE,FALSE),
         v_ep4ew=if_else(ceiling(t_ep4ew)==0&transfer!=0,TRUE,FALSE),
         v_equality=if_else(ceiling(t_equality)==0&transfer!=0,TRUE,FALSE))
```

```{r}
ob_A_lab_long <- ob_A_lab %>%
  mutate(income_cue=if_else(yo_l_01>=yo_h_01 & yo_l_02>=yo_h_02 &yo_l_03>=yo_h_03 &yo_l_04>=yo_h_04 &yo_l_05>=yo_h_05 &yo_l_06>=yo_h_06 &yo_l_07>=yo_h_07 &yo_l_08>=yo_h_08 &yo_l_09>=yo_h_09 &yo_l_10>=yo_h_10,1,0),
         effort_cue=if_else(yo_l_02>=yo_h_01 & yo_l_04>=yo_h_02 & yo_l_06>=yo_h_03 & yo_l_08>=yo_h_04 & yo_l_10>=yo_h_05,1,0),
         increasing = if_else(monotonically_increasing(c(yo_l_00,yo_l_01,yo_l_02,yo_l_03,yo_l_04,yo_l_05,yo_l_06,yo_l_07,yo_l_08,yo_l_09,yo_l_10))&monotonically_increasing(c(yo_h_00,yo_h_01,yo_h_02,yo_h_03,yo_h_04,yo_h_05,yo_h_06,yo_h_07,yo_h_08,yo_h_09,yo_h_10)),1,0),
         decreasing = if_else(monotonically_decreasing(c(yo_l_00,yo_l_01,yo_l_02,yo_l_03,yo_l_04,yo_l_05,yo_l_06,yo_l_07,yo_l_08,yo_l_09,yo_l_10))&monotonically_decreasing(c(yo_h_00,yo_h_01,yo_h_02,yo_h_03,yo_h_04,yo_h_05,yo_h_06,yo_h_07,yo_h_08,yo_h_09,yo_h_10)),1,0),
         constant = if_else(increasing==decreasing & increasing!=0,1,0),
         increasing = increasing - constant,
         decreasing=decreasing-constant,
         all_zero = if_else(constant == 1 & yo_l_00==0,1,0),
         pattern = case_when(increasing==1~"increasing",decreasing==1~"decreasing",all_zero==1~"all_zero",TRUE ~ "other"))%>%  
  pivot_longer(cols=yo_l_00:yo_h_10,names_to = "question",values_to ="transfer") %>%
  mutate(wage_j = if_else(str_sub(question,4,4)=="l",1,2),
         effort_j = as.numeric(str_sub(question,-2)),
         income_j = wage_j*effort_j,
         income_post_transfer = earning-transfer,
         transfer_per=(transfer/earning)*100,
                                belief = case_when(str_sub(question,-2)=="00" & str_sub(question,4,4)=="l" ~believe_income_00,
                                                   str_sub(question,-2)=="01" & str_sub(question,4,4)=="l" ~believe_income_01,
                                                   str_sub(question,-2)=="02" & str_sub(question,4,4)=="l" ~believe_income_02,
                                                   str_sub(question,-2)=="03" & str_sub(question,4,4)=="l" ~believe_income_03,
                                                   str_sub(question,-2)=="04" & str_sub(question,4,4)=="l" ~believe_income_04,
                                                   str_sub(question,-2)=="05" & str_sub(question,4,4)=="l" ~believe_income_05,
                                                   str_sub(question,-2)=="06" & str_sub(question,4,4)=="l" ~believe_income_06,
                                                   str_sub(question,-2)=="07" & str_sub(question,4,4)=="l" ~believe_income_07,
                                                   str_sub(question,-2)=="08" & str_sub(question,4,4)=="l" ~believe_income_08,
                                                   str_sub(question,-2)=="09" & str_sub(question,4,4)=="l" ~believe_income_09,
                                                   str_sub(question,-2)=="10" & str_sub(question,4,4)=="l" ~believe_income_10,
                                                   str_sub(question,-2)=="00" & str_sub(question,4,4)=="h" ~believe_income_00,
                                                   str_sub(question,-2)=="01" & str_sub(question,4,4)=="h" ~believe_income_02,
                                                   str_sub(question,-2)=="02" & str_sub(question,4,4)=="h" ~believe_income_04,
                                                   str_sub(question,-2)=="03" & str_sub(question,4,4)=="h" ~believe_income_06,
                                                   str_sub(question,-2)=="04" & str_sub(question,4,4)=="h" ~believe_income_08,
                                                   str_sub(question,-2)=="05" & str_sub(question,4,4)=="h" ~believe_income_10,
                                                   str_sub(question,-2)=="06" & str_sub(question,4,4)=="h" ~believe_income_12,
                                                   str_sub(question,-2)=="07" & str_sub(question,4,4)=="h" ~believe_income_14,
                                                   str_sub(question,-2)=="08" & str_sub(question,4,4)=="h" ~believe_income_16,
                                                   str_sub(question,-2)=="09" & str_sub(question,4,4)=="h" ~believe_income_18,
                                                   str_sub(question,-2)=="10" & str_sub(question,4,4)=="h" ~believe_income_20,),
         ep4ew=(earning+income_j)/(num_correct+effort_j)*num_correct,
         t_ep4ew=if_else(earning-ep4ew>0,earning-ep4ew,0),
         b_ep4ew=if_else(((floor(t_ep4ew)>0|ceiling(t_ep4ew)>0)&transfer>0)|((floor(t_ep4ew)==0|ceiling(t_ep4ew)==0)&transfer==0),TRUE,FALSE),
         equality=(earning+income_j)/2,
         t_equality=if_else(earning-equality>0,earning-equality,0),
         b_equality=if_else(((floor(t_equality)>0|ceiling(t_equality)>0)&transfer>0)|((floor(t_equality)==0|ceiling(t_equality)==0)&transfer==0),TRUE,FALSE),
         v_ep4ew=if_else(ceiling(t_ep4ew)==0&transfer!=0,TRUE,FALSE),
         v_equality=if_else(ceiling(t_equality)==0&transfer!=0,TRUE,FALSE))
```

## Anomaly
 
Table 1 MTurk
```{r}
ob_A_mTurk_long %>% group_by(code) %>% summarise(n=n(),
                                                v_ep4ew=sum(v_ep4ew)!=0,
                                                v_equality=sum(v_equality)!=0,
                                                v_effort_cue=!effort_cue[1],
                                                v_income_cue=!income_cue[1],
                                                v_deservingness=v_effort_cue|v_income_cue) %>%
  ungroup()%>%
  summarise(nv_ep4ew=sum(v_ep4ew),
            nv_equality=sum(v_equality),
            nv_effort_cue=sum(v_effort_cue),
            nv_income_cue=sum(v_income_cue),
            nv_deservingness=sum(v_deservingness))
```
Table 1 student
```{r}
ob_A_lab_long %>% group_by(label) %>% summarise(n=n(),
                                                v_ep4ew=sum(v_ep4ew)!=0,
                                                v_equality=sum(v_equality)!=0,
                                                v_effort_cue=!effort_cue[1],
                                                v_income_cue=!income_cue[1],
                                                v_deservingness=v_effort_cue|v_income_cue) %>%
  ungroup()%>%
  summarise(nv_ep4ew=sum(v_ep4ew),
            nv_equality=sum(v_equality),
            nv_effort_cue=sum(v_effort_cue),
            nv_income_cue=sum(v_income_cue),
            nv_deservingness=sum(v_deservingness))
```


### MTurk

```{r}
MTurk_Ob_A_iv_intercept <- FLXMRglmfix(~income_j+effort_j,family = "gaussian")
MTurk_Ob_A_ratio_iv_intercept <- stepFlexmix(transfer_per ~ 1 |code, k=1:8, model=MTurk_Ob_A_iv_intercept, data=ob_A_mTurk_long,nrep=30)
```

```{r}
eva_matrix_MTurk_Ob_A_ratio_iv_intercept <- FMR_eva(MTurk_Ob_A_ratio_iv_intercept,ob_A_mTurk_long$code)
eva_matrix_MTurk_Ob_A_ratio_iv_intercept
```

```{r}
MTurk_Ob_A_ratio_iv_intercept_k4 <- refit(MTurk_Ob_A_ratio_iv_intercept@models$`4`,method="mstep")
```


```{r}
MTurk_Ob_A_ratio_iv_intercept_k4_EM <- flexmix(transfer_per ~ 1 |code,
                                         model=Ob_A_iv_intercept,
                                         cluster=posterior(MTurk_Ob_A_ratio_iv_intercept@models$`4`)[ob_A_mTurk_long$pattern!="all_zero",],
                                         data=ob_A_mTurk_long%>%filter(pattern!="all_zero"))
```



```{r}
MTurk_Ob_A_ratio_iv_intercept_k4_EM_refit <- refit(MTurk_Ob_A_ratio_iv_intercept_k4_EM)
summary(MTurk_Ob_A_ratio_iv_intercept_k4_EM_refit)
```

```{r}
table(ob_A_mTurk$pattern,unique(cbind(ob_A_mTurk_long$code,MTurk_Ob_A_ratio_iv_intercept@models$`4`@cluster))[,2])
```

```{r}
table(ob_A_mTurk$type)
```


```{r}
ob_A_mTurk_long$type <- MTurk_Ob_A_ratio_iv_intercept@models$`4`@cluster
```

```{r}
ob_A_mTurk_long <- ob_A_mTurk_long %>% mutate(type=factor(case_when(type==1~"income",type==2~"selfinterest",
                                              type==3~"constant",type==4~"effort"),levels = c("effort","income","selfinterest","constant")))
```

```{r}
ob_A_mTurk_ratio_iv_intercept_regest_type <- data.frame(type=c("effort","income","selfinterest","constant"),
                                                  b_0=NA,
                                                  b_income_j=NA,
                                                  b_effort_j=NA)
ob_A_mTurk_ratio_iv_intercept_regest_type[1,2:4] <- MTurk_Ob_A_ratio_iv_intercept_k4@components[[1]]$Comp.4$coefficients
ob_A_mTurk_ratio_iv_intercept_regest_type[2,2:4] <- MTurk_Ob_A_ratio_iv_intercept_k4@components[[1]]$Comp.1$coefficients
ob_A_mTurk_ratio_iv_intercept_regest_type[3,2:4] <- MTurk_Ob_A_ratio_iv_intercept_k4@components[[1]]$Comp.2$coefficients
ob_A_mTurk_ratio_iv_intercept_regest_type[4,2:4] <- MTurk_Ob_A_ratio_iv_intercept_k4@components[[1]]$Comp.3$coefficients
```

```{r}
ob_A_mTurk_ratio_iv_intercept_regest_type <- ob_A_mTurk_ratio_iv_intercept_regest_type %>%
  mutate(me_wage_low=b_income_j+b_effort_j*10,
         me_wage_high=b_income_j+0.5*b_effort_j*10)
```

```{r}
ob_A_mTurk_long <- ob_A_mTurk_long %>% left_join(ob_A_mTurk_ratio_iv_intercept_regest_type,by="type")
```

### Student

```{r}
Ob_A_iv_intercept <- FLXMRglmfix(~income_j+effort_j,family = "gaussian")
Ob_A_ratio_iv_intercept <- stepFlexmix(transfer_per ~ 1 |label, k=1:8, model=Ob_A_iv_intercept, data=ob_A_lab_long,nrep=30)
```

```{r}
eva_matrix_Ob_A_ratio_iv_intercept <- FMR_eva(Ob_A_ratio_iv_intercept,ob_A_lab_long$label)
```

Table S1 student
```{r}
eva_matrix_Ob_A_ratio_iv_intercept
```


```{r}
Ob_A_ratio_iv_intercept_k4 <- refit(Ob_A_ratio_iv_intercept@models$`4`,method="mstep")
```

```{r}
summary(Ob_A_ratio_iv_intercept_k4)
```

```{r}
Ob_A_ratio_iv_intercept_k4_EM <- flexmix(transfer_per ~ 1 |label,
                                         model=Ob_A_iv_intercept,
                                         cluster=posterior(Ob_A_ratio_iv_intercept@models$`4`)[ob_A_lab_long$pattern!="all_zero",],
                                         data=ob_A_lab_long%>%filter(pattern!="all_zero"))
```

```{r}
table(ob_A_lab$pattern,unique(cbind(ob_A_lab_long$label,Ob_A_ratio_iv_intercept@models$`4`@cluster))[,2])
table(unique(cbind(ob_A_lab_long$label,Ob_A_ratio_iv_intercept@models$`4`@cluster))[,2])
```

```{r}
table(ob_A_lab$type)
```


```{r}
ob_A_lab_long$type <- Ob_A_ratio_iv_intercept@models$`4`@cluster
```

```{r}
ob_A_lab_long <- ob_A_lab_long %>% mutate(type=factor(case_when(type==1~"effort",type==2~"income",
                                              type==3~"selfinterest",type==4~"constant"),levels = c("effort","income","selfinterest","constant")))
```

```{r}
Ob_A_ratio_iv_intercept_regest_type <- data.frame(type=c("effort","income","selfinterest","constant"),
                                                  b_0=NA,
                                                  b_income_j=NA,
                                                  b_effort_j=NA)
Ob_A_ratio_iv_intercept_regest_type[1,2:4] <- Ob_A_ratio_iv_intercept_k4@components[[1]]$Comp.1$coefficients
Ob_A_ratio_iv_intercept_regest_type[2,2:4] <- Ob_A_ratio_iv_intercept_k4@components[[1]]$Comp.2$coefficients
Ob_A_ratio_iv_intercept_regest_type[3,2:4] <- Ob_A_ratio_iv_intercept_k4@components[[1]]$Comp.3$coefficients
Ob_A_ratio_iv_intercept_regest_type[4,2:4] <- Ob_A_ratio_iv_intercept_k4@components[[1]]$Comp.4$coefficients
```

```{r}
Ob_A_ratio_iv_intercept_regest_type <- Ob_A_ratio_iv_intercept_regest_type %>%
  mutate(me_wage_low=b_income_j+b_effort_j,
         me_wage_high=b_income_j+0.5*b_effort_j)
```

### Figures 

Figure 3
```{r}
# Redistribution by types
ob_A_lab_long %>% mutate(sample="Student") %>% 
  bind_rows(ob_A_mTurk_long %>% mutate(sample="MTurk")) %>%
  mutate(wage_level=factor(if_else(wage_j==1 | wage_j==0.1,"low","high")),
         type=fct_relevel(fct_recode(type,`Highly Responsive\nto effort`="effort",`Less responsive\nto effort`="income",`Self-interested`="selfinterest",`Other`="constant"),"Highly Responsive\nto effort","Less responsive\nto effort","Self-interested","Other")) %>%
  ggplot() +
  geom_jitter(aes(x=income_j,y=transfer_per,colour=as.factor(wage_level)),size=0.7)+
  # geom_smooth(aes(x=income_j,y=transfer_per,linetype=as.factor(wage_j),colour=label),method = lm, se = FALSE)+
  geom_abline(aes(slope=me_wage_low,intercept = b_0),colour="#F2C45F",linewidth=1.2,alpha=0.6)+
  geom_abline(aes(slope=me_wage_high,intercept = b_0),colour="#1A80BB",linewidth=1.2,alpha=0.6)+
  scale_x_continuous(name="Recipient's income")+
  scale_y_continuous(name="Percentage of income transferred")+
  scale_color_manual(name="Recipient's wage level",values=c("#1A80BB","#F2C45F"))+
  facet_grid(type~sample, scales = "free_x")+
  theme_bw(base_size=14)+
  theme(legend.position="bottom",strip.text = element_text(size = 14))
# ggsave("Redist_typeXsample.eps",width=88,height=140,units = "mm")
ggsave("./Figures/Figure 3.pdf",width=180,height=225,units = "mm")
```

## Randon-effect Models

### MTurk

```{r}
Ob_A_lme_mTurk <-lmer(transfer_per~income_j+effort_j + (1 | code),
                         data=ob_A_mTurk_long)
summary(Ob_A_lme_mTurk)
```

```{r}
Ob_A_lme_mTurk_no_self <-lmer(transfer_per~income_j+effort_j + (1 | code),
                         data=ob_A_mTurk_long%>%filter(type!="selfinterest"))
summary(Ob_A_lme_mTurk_no_self)
```

### Students

```{r}
Ob_A_lme_lab <-lmer(transfer_per~income_j+effort_j + (1 | label),
                         data=ob_A_lab_long)
summary(Ob_A_lme_lab)
```

```{r}
Ob_A_lme_lab_no_self <-lmer(transfer_per~income_j+effort_j + (1 | label),
                         data=ob_A_lab_long%>%filter(type!="selfinterest"))
summary(Ob_A_lme_lab_no_self)
```

## Multinomial Logit

```{r}
delta.method <- function(func1,beta,V,...)
{
	c_beta <- jacobian(func1,beta,method="simple",...)
	Vc <- c_beta %*% V %*% t(c_beta)
	est <- func1(beta,...)
	se <- sqrt(diag(Vc))
	p <- pnorm(abs(est)/se,lower.tail = FALSE)*2
	list(est=est,se=se,p=p)
}
OR.mlogit <- function(beta){
  exp(beta)
}
cp.mlogit <- function(beta,K,X0,X1){
  beta_mat <- matrix(beta,
                 nrow=K,
                 byrow = TRUE)
  beta_mat <- cbind(matrix(rep(0,K),nrow=K),beta_mat)
  e_mat0 <- exp(X0%*%beta_mat)
  p0 <- t(e_mat0)%*%(1/matrix(rowSums(e_mat0)))/nrow(X0)*100
  e_mat1 <- exp(X1%*%beta_mat)
  p1 <- t(e_mat1)%*%(1/matrix(rowSums(e_mat1)))/nrow(X0)*100
  p1-p0
}
me.mlogit <- function(beta,K,X,k){
  beta_mat <- matrix(beta,
                 nrow=K,
                 byrow = TRUE)
  beta_mat <- cbind(matrix(rep(0,K),nrow=K),beta_mat)
  e_mat <- exp(X%*%beta_mat)
  e_sum <- rowSums(e_mat)
  de <- e_mat*matrix(rep(beta_mat[k,],nrow(X)),nrow=nrow(X),byrow = TRUE)
  de_sum <- rowSums(de)
  colMeans((de*matrix(rep(e_sum,ncol(beta_mat)),nrow=nrow(X))-e_mat*matrix(rep(de_sum,ncol(beta_mat)),nrow=nrow(X)))/matrix(rep(e_sum,ncol(beta_mat)),nrow=nrow(X))^2)*100
}
eff.mlogit <- function(obj,me=c()){
  category <- names(obj$freq)
  L <- length(category)
  X_raw <- data.frame(obj$model)
  K_raw <- ncol(X_raw)-4
  X_raw <- X_raw[X_raw$type==TRUE,1:K_raw+1]
  N <- nrow(X_raw)
  var <- str_unique(str_extract(names(coef(obj)),"^[^:]+"))
  K <- length(var)
  X <- matrix(rep(NA,N*K),nrow=N)
  X[,1] <- 1
  var_ft <- rep(NA,K-1)
  var_grouped <- rep(NA,K-1)
  k <- 1
  
  for(i in 1:K_raw){
    if(is.factor(X_raw[,i])){
      levels <- levels(X_raw[,i])
      var_ft[k:(k+length(levels)-2)] <- TRUE
      var_grouped[k:(k+length(levels)-2)] <- i+1
      for(l in levels[2:length(levels)]){
        k <- k+1
        X[,k] <- X_raw[,i]==l
        
      }
    }else{
      var_ft[k] <- FALSE
      var_grouped[k] <- i+1
      k <- k+1
      X[,k] <- X_raw[,i]
      
    }
  }
  var_ft <- c(FALSE,var_ft)
  var_grouped <- c(1,var_grouped)
  output <- data.frame(var=var)

  for(i in 1:K){
    if(var_ft[i]==FALSE){
      if(var[i] %in% me){
        r <- delta.method(me.mlogit,coef(obj),vcov(obj),K=K,X=X,k=i)
      }else{
        X0 <- X
        X1 <- X
        X1[,i] <- X1[,i]+1
        r <- delta.method(cp.mlogit,coef(obj),vcov(obj),K=K,X0=X0,X1=X1)
      }
    }else{
      X0 <- X
      X0[,var_grouped==var_grouped[i]] <- 0
      X1 <- X0
      X1[,i] <- 1
      r <- delta.method(cp.mlogit,coef(obj),vcov(obj),K=K,X0=X0,X1=X1)
    }
    output[i,2:(1+L)] <- str_c(str_c(round(r$est,3),case_when(r$p<=0.001~"***",
                                                   r$p>0.001&r$p<=0.01~"**",
                                                   r$p>0.01&r$p<=0.05~"*",
                                                   r$p>0.05&r$p<=0.1~"+",
                                                   .default="")),
                                   str_c("(",round(r$se,3),")"),sep = "\n")
  }
  colnames(output)[2:(L+1)] <- category
  output <- output[-1,]
  output<-flextable(output)
  output
}

pp.mlogit <- function(beta,K,X0){
  beta_mat <- matrix(beta,
                 nrow=K,
                 byrow = TRUE)
  beta_mat <- cbind(matrix(rep(0,K),nrow=K),beta_mat)
  e_mat0 <- exp(X0%*%beta_mat)
  # p0 <- e_mat0/sum(e_mat0)
  # p0
  p0 <- t(e_mat0)%*%(1/matrix(rowSums(e_mat0)))/nrow(X0)*100  
  p0
}

X.mlogit <- function(obj){
  category <- names(obj$freq)
  L <- length(category)
  X_raw <- data.frame(obj$model)
  K_raw <- ncol(X_raw)-4
  X_raw <- X_raw[X_raw$type==TRUE,1:K_raw+1]
  N <- nrow(X_raw)
  var <- str_unique(str_extract(names(coef(obj)),"^[^:]+"))
  K <- length(var)
  X <- matrix(rep(NA,N*K),nrow=N)
  X[,1] <- 1
  var_ft <- rep(NA,K-1)
  var_grouped <- rep(NA,K-1)
  k <- 1
  
  for(i in 1:K_raw){
    if(is.factor(X_raw[,i])){
      levels <- levels(X_raw[,i])
      var_ft[k:(k+length(levels)-2)] <- TRUE
      var_grouped[k:(k+length(levels)-2)] <- i+1
      for(l in levels[2:length(levels)]){
        k <- k+1
        X[,k] <- X_raw[,i]==l
        
      }
    }else{
      var_ft[k] <- FALSE
      var_grouped[k] <- i+1
      k <- k+1
      X[,k] <- X_raw[,i]
      
    }
  }
  colnames(X) <- var
  list(X=X,var=var,category=category)
}
```

### MTurk

```{r}
ob_A_mTurk_mnl <- dfidx(ob_A_mTurk %>% 
  mutate(wage_lev=factor(wage,levels=c(0.1,0.2),labels = c("low","high")),
         pct_high=factor(if_else(more_high_wage==1,"75%","25%"),levels=c("25%","75%")),
         college=if_else(education>4,TRUE,FALSE)),
  shape="wide",choice="type")
```

```{r}
ob_A_mTurk_mnp_dm1 <- mlogit(type ~ 0| wage_lev+pct_high+as.factor(gender)+ideo_centered+pid_3point+age+I(log(1+household_income))+college, ob_A_mTurk_mnl,reflevel="effort")
summary(ob_A_mTurk_mnp_dm1)
```

```{r}
flextable(tidy(ob_A_mTurk_mnp_dm1) %>% 
            mutate(OR.p.value=2*pnorm(abs((exp(estimate)-1)/(exp(estimate)*std.error)),lower.tail = FALSE))%>%
  mutate(output=str_c(round(exp(estimate),3),
                            case_when(OR.p.value<=0.001~"***",
                                      OR.p.value<=0.01&p.value>0.001~"**",
                                      OR.p.value<=0.05&p.value>0.01~"*",
                                      OR.p.value<=0.1&p.value>0.05~"+",
                                      OR.p.value>0.1~""),
                      "\n","(",
                      round(exp(estimate)*std.error,3),")"))%>%
  select(term,output))
```

### Students

```{r}
ob_A_lab_mnl <- dfidx(ob_A_lab%>%
  mutate(wage_lev=factor(wage,levels=c(1,2),labels = c("low","high")),
         pct_high=factor(if_else(more_high_wage==1,"75%","25%"),levels=c("25%","75%")),
         hh=if_else(wage_lev=="high"&more_high_wage==1,TRUE,FALSE),
         ideo_none=factor(ideo_none)),shape="wide",choice="type")
```

```{r}
ob_A_lab_mnp_dm1 <- mlogit(type ~ 0| wage_lev+pct_high+as.factor(gender)+ideo_centered+ideo_none+pid_3point+I(log(1+household_income)), ob_A_lab_mnl,reflevel="effort")
summary(ob_A_lab_mnp_dm1)
```

```{r}
flextable(tidy(ob_A_lab_mnp_dm1) %>% 
            mutate(OR.p.value=2*pnorm(abs((exp(estimate)-1)/(exp(estimate)*std.error)),lower.tail = FALSE))%>%
  mutate(output=str_c(round(exp(estimate),3),
                            case_when(OR.p.value<=0.001~"***",
                                      OR.p.value<=0.01&p.value>0.001~"**",
                                      OR.p.value<=0.05&p.value>0.01~"*",
                                      OR.p.value<=0.1&p.value>0.05~"+",
                                      OR.p.value>0.1~""),
                      "\n","(",
                      round(exp(estimate)*std.error,3),")"))%>%
  select(term,output))
```

### Figures

```{r}
s_OR_mnl <- delta.method(OR.mlogit,coef(ob_A_lab_mnp_dm1),vcov(ob_A_lab_mnp_dm1))
m_OR_mnl <- delta.method(OR.mlogit,coef(ob_A_mTurk_mnp_dm1),vcov(ob_A_mTurk_mnp_dm1))
```


```{r}
data.frame(varXtype=names(s_OR_mnl$est),
           OR=s_OR_mnl$est,
           ub=s_OR_mnl$est+qnorm(0.975)*s_OR_mnl$se,
           lb=s_OR_mnl$est+qnorm(0.025)*s_OR_mnl$se,
           sample="Student")%>% rowwise() %>% 
  mutate(var=str_split(varXtype,":")[[1]][1],
         type=str_split(varXtype,":")[[1]][2])%>%
  filter(str_detect(varXtype,"education|Intercept|constant",negate=TRUE))%>%
  ggplot()+
  geom_point(aes(x=type,y=OR))+
  geom_errorbar(aes(x=type,ymin=lb,ymax=ub))+
  geom_hline(aes(yintercept = 1),linetype="dashed")+
  coord_flip()+
  facet_wrap(~var,scales="free_x")
```

Figure S1
```{r}
data.frame(varXtype=names(m_OR_mnl$est),
           OR=m_OR_mnl$est,
           ub=m_OR_mnl$est+qnorm(0.975)*m_OR_mnl$se,
           lb=m_OR_mnl$est+qnorm(0.025)*m_OR_mnl$se,
           sample="MTurk")%>% rowwise() %>% 
  mutate(var=str_split(varXtype,":")[[1]][1],
         type=str_split(varXtype,":")[[1]][2])%>%
  filter(str_detect(varXtype,"education|Intercept|constant",negate=TRUE))%>%
  ggplot()+
  geom_point(aes(x=type,y=OR))+
  geom_errorbar(aes(x=type,ymin=lb,ymax=ub),width=0.5)+
  geom_hline(aes(yintercept = 1),linetype="dashed")+
  coord_flip()+
  facet_wrap(~var,scales="free_x")
```
```{r}
data.frame(varXtype=names(s_OR_mnl$est),
           OR=s_OR_mnl$est,
           ub=s_OR_mnl$est+qnorm(0.975)*s_OR_mnl$se,
           lb=s_OR_mnl$est+qnorm(0.025)*s_OR_mnl$se,
           sample="Student")%>% 
  bind_rows(data.frame(varXtype=names(m_OR_mnl$est),
           OR=m_OR_mnl$est,
           ub=m_OR_mnl$est+qnorm(0.975)*m_OR_mnl$se,
           lb=m_OR_mnl$est+qnorm(0.025)*m_OR_mnl$se,
           sample="MTurk")) %>%
  rowwise() %>% 
  mutate(var=str_split(varXtype,":")[[1]][1],
         type=str_split(varXtype,":")[[1]][2])%>%
  filter(str_detect(varXtype,"Intercept|constant|none",negate=TRUE))%>%
  mutate(var=factor(var,
                    levels = c("wage_levhigh",
                               "pct_high75%",
                               "as.factor(gender)1",
                               "ideo_centered",
                               "pid_3pointInd",
                               "pid_3pointRep",
                               "I(log(1 + household_income))",
                               "collegeTRUE",
                               "age"),
                    labels=c("High wage",
                             "75% high wage",
                             "Male",
                             "Political ideology",
                             "Independent",
                             "Republican",
                             "Log household income",
                             "College degree",
                             "Age")),
         type=factor(type,levels=c("selfinterest","income"),labels=c("Self-interested","Less responsive\nto effort")))%>%
  ggplot()+
  geom_point(aes(x=type,y=OR))+
  geom_errorbar(aes(x=type,ymin=lb,ymax=ub),width=0.3)+
  geom_hline(aes(yintercept = 1),linetype="dashed")+
  scale_x_discrete(name="Type")+
  scale_y_continuous(name="Odds ratio\n(relative to the highly responsive to effort type)")+
  coord_flip()+
  facet_wrap(sample~var,scales="free_x",ncol = 4)+
  theme_bw(base_size=12)+
  theme(strip.text = element_text(size = 10),
        panel.spacing.x = unit(8, "mm"))
ggsave("p_MNL_OR.jpeg",width=240,height=180,units = "mm")
```


```{r}
X <- X.mlogit(ob_A_lab_mnp_dm1)$X
category <- X.mlogit(ob_A_lab_mnp_dm1)$category
pp_lab_by_ideo <- data.frame(ideology=NULL,category=NULL,probability=NULL)
for(i in 0:6){
  X[,"ideo_centered"] <- i
  temp <- pp.mlogit(coef(ob_A_lab_mnp_dm1),ncol(X),X0=X)
  
  pp_temp <- data.frame(ideology=i,category=category,probability=temp)
  pp_lab_by_ideo <- rbind(pp_lab_by_ideo,pp_temp)
}
```

```{r}
X <- X.mlogit(ob_A_mTurk_mnp_dm1)$X
category <- X.mlogit(ob_A_mTurk_mnp_dm1)$category
pp_mturk_by_ideo <- data.frame(ideology=NULL,category=NULL,probability=NULL)
for(i in 0:6){
  X[,"ideo_centered"] <- i
  temp <- pp.mlogit(coef(ob_A_mTurk_mnp_dm1),ncol(X),X0=X)
  
  pp_temp <- data.frame(ideology=i,category=category,probability=temp)
  pp_mturk_by_ideo <- rbind(pp_mturk_by_ideo,pp_temp)
}
```

Figure 4
```{r}
pp_lab_by_ideo %>% mutate(sample="Student") %>%
  bind_rows(pp_mTurk_by_ideo %>% mutate(sample="MTurk"))%>%
  mutate(category=fct_relevel(fct_recode(category,
                                         "Highly responsive\nto effort"="effort",
                                         "Less responsive\nto effort"="income",
                                         "Self-interested"="selfinterest",
                                         "Other"="constant"),
                              "Highly responsive\nto effort",
                              "Less responsive\nto effort",
                              "Self-interested",
                              "Other"))%>%
  ggplot()+
  geom_line(aes(x=ideology,y=probability,colour=category),linewidth=0.7)+
  scale_x_continuous(name="Ideology (from 0 (extremely liberal) to 6 (extremely conservative)")+
  scale_y_continuous(name="Mean predicted probability")+
  facet_wrap(~sample)+
  theme_bw(base_size=14)+
  theme(legend.position="bottom",plot.background = element_rect(fill = "#FAFAFA"),
        legend.title=element_blank())
ggsave("./Figures/Figure 4.pdf",width=180,height=112,units = "mm")

```

## Post Questionnaire

Figure 5
```{r}
ob_A_lab %>%
  mutate(post_q1=(post_q1-1)/9)%>%
  group_by(type)%>%
  summarise(mean=mean(post_q1),
            se=sqrt(var(post_q1)/n()))%>%
  mutate(sample="Student") %>%
  bind_rows(ob_A_mTurk %>%
              mutate(post_q1=(post_q1-1)/9)%>%
  group_by(type)%>%
  summarise(mean=mean(post_q1),
            se=sqrt(var(post_q1)/n()))%>%
  mutate(sample="MTurk")) %>%
  mutate(sample=factor(sample,levels=c("MTurk","Student")),
         type=fct_relevel(fct_recode(type,
                                         "Highly\nresponsive\nto effort"="effort",
                                         "Less\nresponsive\nto effort"="income",
                                         "Self\ninterested"="selfinterest",
                                         "Other"="constant"),
                              "Highly\nresponsive\nto effort",
                              "Less\nresponsive\nto effort",
                              "Self\ninterested",
                              "Other"))%>%
  ggplot(aes(x=type,y=mean,fill=sample))+
  geom_bar(stat="identity",position = "dodge")+
  geom_errorbar(aes(ymin=mean+qnorm(0.025)*se, ymax=mean+qnorm(0.975)*se), width=.2,position=position_dodge(.9)) +
  # geom_text(aes(x=type,y=0.05,label = round(mean,3),fill=sample), colour = "black",size = 4)+
  scale_x_discrete(name="Type")+
  # scale_y_continuous(name="Proportion")+
  scale_fill_manual(name="Sample",values=c("#1A80BB","#F2C45F"))+
  theme_bw(base_size=14)+
  theme(panel.border = element_blank(), 
        panel.spacing.x = unit(0,"line"),
        legend.position="bottom")
ggsave("./Figures/Figure 5.pdf",device="pdf",width=1920,height=1280,units = "px")
```

Figure S2
```{r}
ob_A_lab %>%
  mutate(post_q2=(post_q2-1)/9)%>%
  group_by(type)%>%
  summarise(mean=mean(post_q2),
            se=sqrt(var(post_q2)/n()))%>%
  mutate(sample="Student") %>%
  bind_rows(ob_A_mTurk %>%
              mutate(post_q2=(post_q2-1)/9)%>%
  group_by(type)%>%
  summarise(mean=mean(post_q2),
            se=sqrt(var(post_q2)/n()))%>%
  mutate(sample="MTurk")) %>%
  mutate(sample=factor(sample,levels=c("MTurk","Student")),
         type=fct_relevel(fct_recode(type,
                                         "Highly\nresponsive\nto effort"="effort",
                                         "Less\nresponsive\nto effort"="income",
                                         "Self\ninterested"="selfinterest",
                                         "Other"="constant"),
                              "Highly\nresponsive\nto effort",
                              "Less\nresponsive\nto effort",
                              "Self\ninterested",
                              "Other"))%>%
  ggplot(aes(x=type,y=mean,fill=sample))+
  geom_bar(stat="identity",position = "dodge")+
  geom_errorbar(aes(ymin=mean+qnorm(0.025)*se, ymax=mean+qnorm(0.975)*se), width=.2,position=position_dodge(.9)) +
  # geom_text(aes(x=type,y=0.05,label = round(mean,3),fill=sample), colour = "black",size = 4)+
  scale_x_discrete(name="Type")+
  # scale_y_continuous(name="Proportion")+
  scale_fill_manual(name="Sample",values=c("#1A80BB","#F2C45F"))+
  theme_bw(base_size=14)+
  theme(panel.border = element_blank(), 
        panel.spacing.x = unit(0,"line"),
        legend.position="bottom")
# ggsave("./Figures/PostQ2_byType.jpeg",device="jpeg",width=1920,height=1280,units = "px")
```

# Income level information only

MTurk
```{r}
no_A_mturk <- read_csv("Unobservable_mturk.csv")
```

Transform to the long format
```{r}
no_A_mturk_long <- no_A_mturk %>% pivot_longer(cols=no_00:no_20,names_to = "question",values_to ="transfer") %>%
  mutate(income_j = round(as.numeric(str_sub(question,-2))*0.1,1),
         wage_j = factor(case_when(income_j %in% c(0.1,0.3,0.5,0.7,0.9)~"low",
                              income_j %in% c(1.2,1.4,1.6,1.8,2)~"high",
                              income_j %in% c(0,0.2,0.4,0.6,0.8,1)~"unknown"),levels=c("unknown","low","high")),
         transfer_per=(transfer/earning)*100,
                                  remain = earning-transfer,
                                  i_larger_income = if_else(earning>income_j,1,0),
                                  belief = case_when(str_sub(question,-2)=="00"~believe_income_00,
                                                     str_sub(question,-2)=="01"~believe_income_01,
                                                     str_sub(question,-2)=="02"~believe_income_02,
                                                     str_sub(question,-2)=="03"~believe_income_03,
                                                     str_sub(question,-2)=="04"~believe_income_04,
                                                     str_sub(question,-2)=="05"~believe_income_05,
                                                     str_sub(question,-2)=="06"~believe_income_06,
                                                     str_sub(question,-2)=="07"~believe_income_07,
                                                     str_sub(question,-2)=="08"~believe_income_08,
                                                     str_sub(question,-2)=="09"~believe_income_09,
                                                     str_sub(question,-2)=="10"~believe_income_10,
                                                     str_sub(question,-2)=="12"~believe_income_12,
                                                     str_sub(question,-2)=="14"~believe_income_14,
                                                     str_sub(question,-2)=="16"~believe_income_16,
                                                     str_sub(question,-2)=="18"~believe_income_18,
                                                     str_sub(question,-2)=="20"~believe_income_20),
         effort_j_ind=case_when(wage_j=="low"~income_j*10,
                                wage_j=="high"~income_j*10/2,
                                wage_j=="unknown"~more_high_wage*(0.75*income_j*10/2+0.25*income_j*10)+
                                  !more_high_wage*(0.25*income_j*10/2+0.75*income_j*10)),
         effort_j_b=belief*0.01*income_j*10/2+(100-belief)*0.01*income_j*10)

```

Student sample
```{r}
no_A_lab <- read_csv("Unobservable_student.csv")
```

```{r}
no_A_lab_long <- no_A_lab %>% pivot_longer(cols=no_00:no_20,names_to = "question",values_to ="transfer") %>%
  mutate(income_j = as.numeric(str_sub(question,-2)),
         wage_j = factor(case_when(income_j %in% c(1,3,5,7,9)~"low",
                              income_j %in% c(12,14,16,18,20)~"high",
                              income_j %in% c(0,2,4,6,8,10)~"unknown"),levels=c("unknown","low","high")),
         transfer_per=(transfer/earning)*100,
                                  remain = earning-transfer,
                                  i_larger_income = if_else(earning>income_j,1,0),
                                  belief = case_when(str_sub(question,-2)=="00"~believe_income_00,
                                                     str_sub(question,-2)=="01"~believe_income_01,
                                                     str_sub(question,-2)=="02"~believe_income_02,
                                                     str_sub(question,-2)=="03"~believe_income_03,
                                                     str_sub(question,-2)=="04"~believe_income_04,
                                                     str_sub(question,-2)=="05"~believe_income_05,
                                                     str_sub(question,-2)=="06"~believe_income_06,
                                                     str_sub(question,-2)=="07"~believe_income_07,
                                                     str_sub(question,-2)=="08"~believe_income_08,
                                                     str_sub(question,-2)=="09"~believe_income_09,
                                                     str_sub(question,-2)=="10"~believe_income_10,
                                                     str_sub(question,-2)=="12"~believe_income_12,
                                                     str_sub(question,-2)=="14"~believe_income_14,
                                                     str_sub(question,-2)=="16"~believe_income_16,
                                                     str_sub(question,-2)=="18"~believe_income_18,
                                                     str_sub(question,-2)=="20"~believe_income_20),
         effort_j_ind=case_when(wage_j=="low"~income_j,
                                wage_j=="high"~income_j/2,
                                wage_j=="unknown"~more_high_wage*(0.75*income_j/2+0.25*income_j)+
                                  !more_high_wage*(0.25*income_j/2+0.75*income_j)),
         effort_j_b=belief*0.01*income_j/2+(100-belief)*0.01*income_j)

```

## Finite mixture model
### MTurk
```{r}
MTurk_No_A_iv_wage_j <- FLXMRglmfix(~as.factor(wage_j)*income_j,family = "gaussian")
MTurk_No_A_ratio_iv_wage_j <- stepFlexmix(transfer_per ~ 1 |code, k=1:8, model=MTurk_No_A_iv_wage_j, data=no_A_mturk_long,nrep=30)
```

```{r}
eva_matrix_MTurk_No_A_ratio_iv_wage_j <- FMR_eva(MTurk_No_A_ratio_iv_wage_j,no_A_mturk_long$code)
eva_matrix_MTurk_No_A_ratio_iv_wage_j
```

```{r}
MTurk_No_A_ratio_iv_wage_j_k5 <- refit(MTurk_No_A_ratio_iv_wage_j@models$`5`,method="mstep")
summary(MTurk_No_A_ratio_iv_wage_j_k5)
```

```{r}
table(unique(cbind(no_A_mturk_long$code,MTurk_No_A_ratio_iv_wage_j@models$`5`@cluster))[,2])
```


```{r}
no_A_mturk_long$type <- MTurk_No_A_ratio_iv_wage_j@models$`5`@cluster
```


```{r}
MTurk_No_A_ratio_iv_wage_j_k5_EM <- flexmix(transfer_per ~ 1 |code,
                                         model=MTurk_No_A_iv_wage_j,
                                         cluster=posterior(MTurk_No_A_ratio_iv_wage_j@models$`5`)[no_A_mturk_long$type!=2,],
                                         data=no_A_mturk_long%>%filter(type!=2))
```

```{r}
MTurk_No_A_ratio_iv_wage_j_k5_EM_refit <- refit(MTurk_No_A_ratio_iv_wage_j_k5_EM)
summary(MTurk_No_A_ratio_iv_wage_j_k5_EM_refit)
```


```{r}
no_A_mturk_long <- no_A_mturk_long %>%
  mutate(type=dplyr::recode(type, `1`="other",`2`="selfinterest",`3`="effort",`4`="income_high",`5`="income_low"))
```



```{r}
MTurk_No_A_ratio_iv_wage_j_k5_EM_refit@components[[1]]$Comp.1[1:6,1]
```


```{r}
MTurk_No_A_ratio_iv_wage_j_regest_type <- data.frame(type=c("other","selfinterest","effort","income_high","income_low"),
                                                  b0=NA,
                                                  b0_l=NA,
                                                  b0_h=NA,
                                                  b_income_j=NA,
                                                  b_income_j_l=NA,
                                                  b_income_j_h=NA)
MTurk_No_A_ratio_iv_wage_j_regest_type[1,2:7] <- MTurk_No_A_ratio_iv_wage_j_k5_EM_refit@components[[1]]$Comp.1[1:6,1]
MTurk_No_A_ratio_iv_wage_j_regest_type[2,2:7] <- 0
MTurk_No_A_ratio_iv_wage_j_regest_type[3,2:7] <- MTurk_No_A_ratio_iv_wage_j_k5_EM_refit@components[[1]]$Comp.2[1:6,1]
MTurk_No_A_ratio_iv_wage_j_regest_type[4,2:7] <- MTurk_No_A_ratio_iv_wage_j_k5_EM_refit@components[[1]]$Comp.3[1:6,1]
MTurk_No_A_ratio_iv_wage_j_regest_type[5,2:7] <- MTurk_No_A_ratio_iv_wage_j_k5_EM_refit@components[[1]]$Comp.4[1:6,1]
```

```{r}
MTurk_No_A_ratio_iv_wage_j_regest_type
```

```{r}
no_A_mturk_long <- no_A_mturk_long %>% left_join(MTurk_No_A_ratio_iv_wage_j_regest_type,by="type")
```

```{r}
no_A_mturk_long <- no_A_mturk_long %>% 
  mutate(e_transfer=b0+b0_l*if_else(wage_j=="low",1,0)+b0_h*if_else(wage_j=="high",1,0)+
           (b_income_j+b_income_j_l*if_else(wage_j=="low",1,0)+b_income_j_h*if_else(wage_j=="high",1,0))*income_j)
```

### Student
```{r}
Lab_No_A_iv_wage_j <- FLXMRglmfix(~as.factor(wage_j)*income_j,family = "gaussian")

Lab_No_A_ratio_iv_wage_j <- stepFlexmix(transfer_per ~ 1 |label, k=1:8, model=Lab_No_A_iv_wage_j, data=no_A_lab_long,nrep=30)
```

```{r}
eva_matrix_Lab_No_A_ratio_iv_wage_j <- FMR_eva(Lab_No_A_ratio_iv_wage_j,no_A_lab_long$label)
eva_matrix_Lab_No_A_ratio_iv_wage_j
```

```{r}
Lab_No_A_ratio_iv_wage_j_k4 <- refit(Lab_No_A_ratio_iv_wage_j@models$`4`,method="mstep")
summary(Lab_No_A_ratio_iv_wage_j_k4)
```


```{r}
table(unique(cbind(no_A_lab_long$label,Lab_No_A_ratio_iv_wage_j@models$`4`@cluster))[,2])
```

```{r}
no_A_lab_long$type <- Lab_No_A_ratio_iv_wage_j@models$`4`@cluster
```

```{r}
Lab_No_A_ratio_iv_wage_j_k4_EM <- flexmix(transfer_per ~ 1 |label,
                                         model=Lab_No_A_iv_wage_j,
                                         cluster=posterior(Lab_No_A_ratio_iv_wage_j@models$`4`)[no_A_lab_long$type!=3,],
                                         data=no_A_lab_long%>%filter(type!=3))
```

```{r}
Lab_No_A_ratio_iv_wage_j_k4_EM_refit <- refit(Lab_No_A_ratio_iv_wage_j_k4_EM)
summary(Lab_No_A_ratio_iv_wage_j_k4_EM_refit)
```

```{r}
no_A_lab_long <- no_A_lab_long %>%
  mutate(type=dplyr::recode(type,`1`="income_high",`2`="income_low",`3`="selfinterest",`4`="other"))
```


```{r}
Lab_No_A_ratio_iv_wage_j_regest_type <- data.frame(type=c("income_high","income_low","selfinterest","other"),
                                                  b0=NA,
                                                  b0_l=NA,
                                                  b0_h=NA,
                                                  b_income_j=NA,
                                                  b_income_j_l=NA,
                                                  b_income_j_h=NA)
Lab_No_A_ratio_iv_wage_j_regest_type[1,2:7] <- Lab_No_A_ratio_iv_wage_j_k4_EM_refit@components[[1]]$Comp.1[1:6,1]
Lab_No_A_ratio_iv_wage_j_regest_type[2,2:7] <-Lab_No_A_ratio_iv_wage_j_k4_EM_refit@components[[1]]$Comp.2[1:6,1]
Lab_No_A_ratio_iv_wage_j_regest_type[3,2:7] <- 0
Lab_No_A_ratio_iv_wage_j_regest_type[4,2:7] <- Lab_No_A_ratio_iv_wage_j_k4_EM_refit@components[[1]]$Comp.3[1:6,1]
```

```{r}
no_A_lab_long <- no_A_lab_long %>% left_join(Lab_No_A_ratio_iv_wage_j_regest_type,by="type")
```

```{r}
no_A_lab_long <- no_A_lab_long %>% 
  mutate(e_transfer=b0+b0_l*if_else(wage_j=="low",1,0)+b0_h*if_else(wage_j=="high",1,0)+
           (b_income_j+b_income_j_l*if_else(wage_j=="low",1,0)+b_income_j_h*if_else(wage_j=="high",1,0))*income_j)
```

### Figures
```{r}
no_A_lab_long %>% group_by(income_j,wage_j,type) %>%
  summarise(mean=mean(transfer_per),
            expected=mean(e_transfer))%>%
  mutate(sample="Student") %>% 
  bind_rows(no_A_mturk_long %>% group_by(income_j,wage_j,type) %>%
  summarise(mean=mean(transfer_per),
            expected=mean(e_transfer))%>%
  mutate(sample="MTurk"))%>% ungroup()%>%
  mutate(type=fct_relevel(fct_recode(type,`Highly responsive\nto income`="income_high",`Less responsive\nto income`="income_low",`Responsive\nto effort`="effort",`Self-interested`="selfinterest",`Other`="other"),"Highly responsive\nto income","Less responsive\nto income","Responsive\nto effort","Self-interested","Other")) %>%
  ggplot()+
  geom_bar(aes(x=as.factor(income_j),y=mean,fill=wage_j),stat="identity",alpha=0.7)+
  # geom_errorbar(aes(x=as.factor(income_j),ymin=lb,ymax=ub))+
  geom_line(aes(x=as.factor(income_j),y=expected,group=wage_j,colour=wage_j),linewidth=1.2)+
  scale_x_discrete(name="Recipient's income")+
  scale_y_continuous(name="Percentage of income transferred")+
  scale_fill_manual(name="Recipient's wage level",values=c("#66c2a5","#F2C45F","#1A80BB"))+
  scale_color_manual(name="Recipient's wage level",values=c("#66c2a5","#F2C45F","#1A80BB"))+
  facet_wrap(sample~type,scales = "free_x",ncol=2,dir = "v")+
  theme_bw(base_size=12)+
  theme(legend.position="bottom",strip.text = element_text(size = 10))
ggsave("./Figures/Figure 6.pdf",width=190,height=270,units = "mm")
```

## Random effect
### MTurk

```{r}
No_A_lme_mTurk <-lmer(transfer_per~as.factor(wage_j)*income_j + (1 | code),
                         data=no_A_mturk_long)
summary(No_A_lme_mTurk)
```

```{r}
No_A_lme_mTurk_no_self <-lmer(transfer_per~as.factor(wage_j)*income_j + (1 | code),
                         data=no_A_mturk_long%>%filter(type!="selfinterest"))
summary(No_A_lme_mTurk_no_self)
```

### Students

```{r}
No_A_lme_lab <-lmer(transfer_per~as.factor(wage_j)*income_j + (1 | label),
                         data=no_A_lab_long)
summary(No_A_lme_lab)
```

```{r}
No_A_lme_lab_no_self <-lmer(transfer_per~as.factor(wage_j)*income_j + (1 | label),
                         data=no_A_lab_long%>%filter(type!="selfinterest"))
summary(No_A_lme_lab_no_self)
```

## Beliefs

Figure S3
```{r}
no_A_mturk_long %>% 
  group_by(more_high_wage,type,wage_j,income_j)%>%
  summarise(mean=mean(belief))%>% ungroup()%>%
  mutate(type=fct_relevel(fct_recode(type,`Highly responsive\nto income`="income_high",`Less responsive\nto income`="income_low",`Responsive\nto effort`="effort",`Self-interested`="selfinterest",`Other`="other"),"Highly responsive\nto income","Less responsive\nto income","Responsive\nto effort","Self-interested","Other"),
         more_high_wage=factor(more_high_wage,labels = c("25% high wage rate", "75% high wage rate"))) %>%
  ggplot()+
  geom_bar(aes(x=as.factor(income_j),y=mean,fill=wage_j),stat="identity")+
  scale_x_discrete(name="Recipient's income")+
  scale_y_continuous(name="Belief of the percentages of high wage rate participants")+
  scale_fill_manual(name="Recipient's wage level",values=c("#66c2a5","#F2C45F","#1A80BB"))+
  facet_grid(type~more_high_wage)+
  theme_bw(base_size=12)+
  theme(legend.position="bottom",strip.text = element_text(size = 10))
ggsave("p_no_MTurk_Belief_typeXsample.jpeg",width=190,height=270,units = "mm")
```

Figure S4
```{r}
no_A_lab_long %>% 
  group_by(more_high_wage,type,wage_j,income_j)%>%
  summarise(mean=mean(belief))%>% ungroup()%>%
  mutate(type=fct_relevel(fct_recode(type,`Highly responsive\nto income`="income_high",`Less responsive\nto income`="income_low",`Self-interested`="selfinterest",`Other`="other"),"Highly responsive\nto income","Less responsive\nto income","Self-interested","Other"),
         more_high_wage=factor(more_high_wage,labels = c("25% high wage rate", "75% high wage rate"))) %>%
  ggplot()+
  geom_bar(aes(x=as.factor(income_j),y=mean,fill=wage_j),stat="identity")+
  scale_x_discrete(name="Recipient's income")+
  scale_y_continuous(name="Belief of the percentages of high wage rate participants")+
  scale_fill_manual(name="Recipient's wage level",values=c("#66c2a5","#F2C45F","#1A80BB"))+
  facet_grid(type~more_high_wage)+
  theme_bw(base_size=12)+
  theme(legend.position="bottom",strip.text = element_text(size = 10))
ggsave("p_no_Student_Belief_typeXsample.jpeg",width=190,height=270,units = "mm")
```

```{r}
no_A_mturk_lm_belief_byMore_high_wage <-lm(belief~ more_high_wage*income_j,
                                 data=no_A_mturk_long%>%
                                   filter(wage_j=="unknown"))
coeftest(no_A_mturk_lm_belief_byMore_high_wage,vcov. = vcovCL(no_A_mturk_lm_belief_byMore_high_wage, cluster=~code))
```

```{r}
no_A_lab_lm_belief_byType <-lm(belief~ more_high_wage*income_j*type,
                                 data=no_A_lab_long%>%
                                   mutate(type=factor(type,levels=c("income_high","income_low","selfinterest","other")))%>%
                                   filter(wage_j=="unknown"))
coeftest(no_A_lab_lm_belief_byType,vcov. = vcovCL(no_A_lab_lm_belief_byType, cluster=~label))
```

Table S10
```{r}
flextable(tidy(no_A_mturk_lm_belief_byType)%>%mutate(MTurk=str_c(round(estimate,3),
                                                        case_when(p.value<=0.001~"***",
                                                                  p.value<=0.01&p.value>0.001~"**",
                                                                  p.value<=0.05&p.value>0.01~"*",
                                                                  p.value<=0.1&p.value>0.05~"+",
                                                                  p.value>0.1~""),
                                                        "\n(",round(std.error,3),")"))%>%
  select(term,MTurk)%>%
  full_join(tidy(no_A_lab_lm_belief_byType)%>%mutate(Student=str_c(round(estimate,3),
                                                        case_when(p.value<=0.001~"***",
                                                                  p.value<=0.01&p.value>0.001~"**",
                                                                  p.value<=0.05&p.value>0.01~"*",
                                                                  p.value<=0.1&p.value>0.05~"+",
                                                                  p.value>0.1~""),
                                                        "\n(",round(std.error,3),")"))%>%
  select(term,Student)))
```

```{r}
no_A_mturk_lm_effort_j_b <-lm(transfer_per~ type*effort_j_b+type*income_j, data=no_A_mturk_long%>%filter(wage_j=="unknown")%>%
                                mutate(type=factor(type,levels=c("income_high","income_low","effort","selfinterest","other"))))
# summary(no_A_mturk_lm_effort_j_b_income_low)
coeftest(no_A_mturk_lm_effort_j_b,vcov. = vcovCL(no_A_mturk_lm_effort_j_b,cluster=~code))
```

```{r}
no_A_lab_lm_effort_j_b <-lm(transfer_per~ type*effort_j_b+type*income_j, data=no_A_lab_long%>%filter(wage_j=="unknown")%>%
                                mutate(type=factor(type,levels=c("income_high","income_low","selfinterest","other"))))
# summary(no_A_mturk_lm_effort_j_b_income_low)
coeftest(no_A_lab_lm_effort_j_b,vcov. = vcovCL(no_A_lab_lm_effort_j_b,cluster=~label))
```

Table S11
```{r}
flextable(tidy(no_A_mturk_lm_effort_j_b)%>%mutate(MTurk=str_c(round(estimate,3),
                                                        case_when(p.value<=0.001~"***",
                                                                  p.value<=0.01&p.value>0.001~"**",
                                                                  p.value<=0.05&p.value>0.01~"*",
                                                                  p.value<=0.1&p.value>0.05~"+",
                                                                  p.value>0.1~""),
                                                        "\n(",round(std.error,3),")"))%>%
  select(term,MTurk)%>%
  full_join(tidy(no_A_lab_lm_effort_j_b)%>%mutate(Student=str_c(round(estimate,3),
                                                        case_when(p.value<=0.001~"***",
                                                                  p.value<=0.01&p.value>0.001~"**",
                                                                  p.value<=0.05&p.value>0.01~"*",
                                                                  p.value<=0.1&p.value>0.05~"+",
                                                                  p.value>0.1~""),
                                                        "\n(",round(std.error,3),")"))%>%
  select(term,Student)))
```

